capture: program drop genatmediant
program genatmediant, rclass
* run after estimation with logistic; generates predicted probability with all variables except varname at median
* By default, restricts to e(sample); this may be overridden by if
* option square: also uses varname2
* option atlist: specifies other variables to set 
* option all: predicts probability for every observation in the sample; by default predicts only once for each unique value   
* REDUNDANT option percentile: also generates binary variable coded 1 at 10th, 50th, 90th percentile 
* REDUNDANT option: sort


syntax varname [if], GENerate(name) [percentile sort SQuare atlist(string) all]
local indeps: colnames e(b)
tempname 	i xb se ci p10 p50 p90 // offset
tempvar		sample

	scalar `ci' = invnormal(.975)
	if "`if'"=="" {
		gen byte `sample' = e(sample) 
		}
	else {
		gen byte `sample' = 0
		replace `sample' = 1 `if'
		local noesample "noesample"
		}
	sort `sample' `varlist' // stack observations not in e(sample) at beginning  
	capture: drop `generate'_xb
	capture: drop `generate'_p
	capture: drop `generate'_ul
	capture: drop `generate'_ll
	quietly gen `generate'_xb=.
	quietly gen `generate'_p=.
	quietly gen `generate'_ul=.
	quietly gen `generate'_ll=.

	quietly margins `if', predict(xb) at(   (median) _all `atlist') `noesample'
	matrix x = r(b)
	local pmedian = 1/(1+exp(-	(x[1,1])	))
	di "Overall probability: " (`pmedian')	
	
	
*	scalar `offset' = 0 // ln(`exposure')
	quietly count if `sample'==0
	di "Observations in sample: " ( _N - r(N) )
	scalar `i' = r(N)+1
	while `i' <= _N {
*		scalar `v' = `varlist'[`i']
*		di `varlist'[`i']
*		local _value = "`v'"
*		di `_value'
*		di `_value'*10
*		di "`varlist'"
		local at = "`varlist'=(" + string(`varlist'[`i']) + ")"
		if "`square'" != "" {
			local at = "`at'" + "  `varlist'2=(" + string(`varlist'2[`i']) + ")"
			}
*		di "`at'" 
		if `i' == 1 | `varlist'[`i']!=`varlist'[`i'-1] {
			quietly margins `if', predict(xb) at(   (median) _all  `atlist' `at'   ) `noesample'
			matrix x = r(b)
			scalar `xb' = x[1,1] //+ `offset'
			quietly replace `generate'_xb = `xb' if _n == `i'
			quietly replace `generate'_p = 1/(1+exp(-	(`xb')	)) if _n == `i'  
			matrix x = r(V)
			scalar `se' = sqrt(x[1,1]) * `ci'
*			di `se'
			quietly replace `generate'_ul = 1/(1+exp(-	(`xb'+`se')	)) if _n == `i'
			quietly replace `generate'_ll = 1/(1+exp(-	(`xb'-`se')	)) if _n == `i'
			di `i' " " _c 
* 			di "  Offset " (`xb'-(`generate'[`i'])) "   Exposure " exp(  (`xb'-(`generate'[`i']))  ) 	
			}
		else if "`all'" != "" {
			quietly replace `generate'_xb = `generate'_xb[`i'-1] if _n == `i'
			quietly replace `generate'_p = 	`generate'_p[`i'-1] if _n == `i'  
			quietly replace `generate'_ul = `generate'_ul[`i'-1] if _n == `i'
			quietly replace `generate'_ll = `generate'_ll[`i'-1] if _n == `i'			
			}
		scalar `i' = `i'+1
		}
*/
	
	if "`percentile'" != "" {
		capture: drop `generate'_percentile
		sort `varlist'
		count 
		scalar `p10' = round(r(N)/10)
		scalar `p50' = round(r(N)/2)
		scalar `p90' = round(r(N)/10*9)
		gen byte `generate'_percentile = 0
		replace `generate'_percentile = 1 if _n==`p10' 	
		replace `generate'_percentile = 1 if _n==`p50' 	
		replace `generate'_percentile = 1 if _n==`p90' 	
		}
		
	capture: drop `generate'_weight
	by `varlist', sort: egen `generate'_weight = total(`sample') 
	
	return scalar pmedian = `pmedian'
end

